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ABSTRACT 

Stellar population synthesis can be approached as an inverse problem. The physical 
information is extracted from the observations through an inverse model. The process 
requires the transformation of the observational errors into model errors. A description 
is given for the error analysis to obtain objectively the errors in the model. 

Finding a solution for overdetermined and under-determined case was the purpose 
of two preceding papers. This new one completes the problem of stellar populations 
synthesis by means of a data base, by providing practical formulae defining the set of 
acceptable solutions. All solutions within this set are compatible, at a given confidence 
level, with the observations. 
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1 INTRODUCTION 

Important astrophysical issues could be solved if we were 
able to deduce the stellar population from the integrated 
light received from a far-away galaxy (or a region of a 
galaxy) . Two methods try to solve this problem : the stellar 
data base synthesis and the evolutionary population synthe- 
sis. The first method is a more empirical approach, it relies 
strictly on the spectral data base of stars or stellar aggre- 
gates, while the second one is a more theoretical approach; it 
relies on our best knowledge on the formation and evolution 
of stars and stellar aggregates. Stellar data base synthesis 
is used in both inverse and forward modeling, while evolu- 
tionary population synthesis is mostly applied as a forward 
model. In our opinion the two methods are complementary 
and should be used jointly in an effort to solve this difficult 
problem. Both methods are subject to errors and ideally 
their results should always be presented with a set of ac- 
ceptable solutions rather than just a solution. The purpose 
of this paper is to establish this set of acceptable solutions 
for the stellar data base synthesis. 

Pelat (1997 and 1998; hereafter paper I and paper II) 
re-investigated the stellar data base synthesis method and 
proposed new and fast algorithms to search for physical so- 
lutions. The sensitivity of this solution to the observational 
errors was discussed in general terms: Monte-Carlo simula- 
tions were performed but no real error analysis had been 
done. This paper provides the user with practical formulae 
which give an estimate of the domain of acceptable solutions 
around the solution. 



In the following § we recall briefly the solution we 
have proposed and then we focus in great detail on the ques- 
tion of error analysis. 



2 POPULATION SYNTHESIS AS AN 
INVERSE PROBLEM 

First we specify what is meant by observables and parame- 
ters in the description of the physical model. 

Observables. It is usually agreed to adopt equivalent 
widths as spectral observables, in order to avoid data re- 
duction errors and to minimize extinction problems. These 
observables are noted here as the vector Wobs. The aim of 
the stellar population synthesis is to find the combinations 
of stars which best reproduce the equivalent widths Wobs of 
the observed object (say, a galaxy). 

Model parameters. The galaxy is described by the un- 
knowns ki. They give the proportional contribution to the 
luminosity due to stars of type i at a reference wavelength 
Ao. The assembly of fc^s is described by the vector k. 

The equation governing the stellar population synthe- 
sis problem is a non linear one, where a set of synthetic 
equivalent widths is a function of a set of stellar luminosity 
contributions ki. We have: 

W.ynj = ^n, , , , J = 1 , . . . , TlA . (1) 

Where, nx is the number of observed equivalent widths; n* 
is the number of stars used for the synthesis; Wji and Iji 
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are respectively the equivalent width and continuum flux at 
wavelength Xj of stars of kind i. The continuum fluxes are 
normalized to one at the reference wavelength Aq. For short 
we write Eq. (0): Wsyn = fik). 

In addition we have the condition that the contribution 
of a star to the luminosity cannot be negative and that the 
normalized proportions ki add up to one. A physical solution 
of Eq.(0) therefore lies in the vector set: 

S = {(fci, . . . , fc„ J I fc. > 0, = !}• (2) 

This set <S is a simplex (e.g. an equilateral triangle for a data 
base of three stars, a tetrahedron for a data base of four and 
so on). The set of all equivalent widths which are able to be 
exactly synthesized by the data base is the image of 5 by ip. 
This is the synthetic domain W introduced in paper I: 

W = {W \ W ^ Lp{k),\/k e S}. (3) 

If Wobs £ W there is at least one exact solution k such that 
ifi{k) = Wobs. 

The set IC of all exact solutions is given formally by: 
K. = v5~^(Wobs). This set may be empty. The solution set 
IC is non-empty if there exists at least one vector k solution 
of the following system: 

fc > 0, ^ 

Afc = 0, I (4) 
b^fc = 1. J 

The matrix elements of A are : [A]^.^ — {Wobsj — Wjijlji and 
is a line of 'ones'. This system may or may not possess a 
solution. 

2.1 Test for the existence of an exact solution 

It is very easy to check if the above system possesses 
at least one solution. If we add to it a linear equation to 
be maximized (e.g. c^k — max!), this system turns into a 
linear program. In fact, it is not necessary to explicit that 
supplementary linear equation, we only have to solve the 
linear program for a feasible solution (see for example Press 
et al. 1992 section 10.8). If no feasible solution exists this 
means that the galaxy cannot be exactly synthesized by the 
data base. This test is extremely rapid. 



2.2 Under-determined case: riA < n* — 1 

If the galaxy can be exactly synthesized (/C is non-empty), 
it has been demonstrated in paper II that /C is a polytope, 
i.e. it is the convex hull of some flnite extreme solutions. A 
solution k is extreme if it possesses at least no — 1 compo- 
nents (stellar contributions) equal zero. The key parameter 
no is the dimension of the null space associated with matrix 
A. The data base is non-degenerate for the observation if we 
have no = n* — nA . 

If the galaxy cannot be synthesized, an approximate 
solution is found with a data base reduced to at most n\ 
stars. The procedure is described in the next section. 



2.3 Overdetermined case: n\ > rii. 

In the overdetermined case, there is most probably no exact 
solution of (|^, the observation is not on the synthetic sur- 
face and one must content oneself with an approximate solu- 
tion. This approximate solution k is usually accepted in the 
least square meaning, that is k must minimize a quadratic 
form : 

= (Wobs - <^(fc))^S-^(Wobs - ^(fe)) , (5) 

where is a positive deflnite 'weight' matrix. It has been 
shown in paper I that a very successful estimate of k is found 
near the 'first-guess' feo: 

feo = (B^B)-'6[b'^(B'^B)-'f)]-' , 

where B is equal to the matrix A augmented with the line 
. It is shown in appendix ^ how one can get rid of the 
positivity constraint fc > when searching for a minimum 
in the neighbourhood of fco. The estimate feo is the unique 
solution of the problem if the matrix is diagonal with 
diagonal elements equal to I^ynj (-^ayn j = X/T-^i ^ji^i)- "^^^ 
argued in paper I that a reasonable weight matrix should not 
be very different from this particular diagonal matrix. 



3 PURPOSES OF THE ERROR ANALYSIS 

An observation of a galaxy is not just the point Wobs in W- 
space, it includes all points interior to an error zone around 
the observation. If we assume the observational errors to be 
Gaussian, this error zone is an hyper-ellipsoid around Wobs- 
We have observed Wobs but we consider that we could have 
observed, with probability 7, any other point at the interior 
of this hyper-ellipsoid. This set of 'probable' points: is 
formally defined as 

£^ = {W\{W~ Wobs)^V-'(W - Wobs) < F-^^-y)} , (6) 

where V is the variance-covariance matrix of the observation 
and F^2 is the distribution function of a variate with n^ 
degrees of freedom. For short, we call S-y the 'error ellipsoid'. 
We recall that V is equal to (AWAW^), where ( ) stands 
for the expectation and AW for a variation of W around 
its mean. 

The purpose of the error analysis is to determine how 
the observational errors are transformed by the inversion 
process. Usually one tries to solve two problems as follows: 
i) deduce Vfc the variance-covariance matrix of the solution 
knowing the matrix V for the observational errors; ii) con- 
struct /C-y, the set of all solutions that fall within the 'error 
ellipsoid', i.e. ICj = ip~^{£^). 

These problems are simplified if the error ellipsoid cor- 
responds to small deviation dW around the observation. 
In that approximation one can linearize (^~^, which implies 
that K.^ is also an hyper-ellipsoid. We have: 

/C, ^ {fe' I (fe' - fc)^P(fe' - fe) < F-/(7)} . (7) 

The matrix P would be equal to VjT^ if were invert- 
ible but this is not the case here: the matrix Vfe is singular 
because of the normalization constraint on fc. Using the def- 
inition of KL-i (^) , one is able to test if an alternative solution 
fc' is acceptable at the confidence level 7. 
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The results of the error analysis are given in 
that, we give in 



Before 



3.1 an overview of the main processes at 
the origin of the observational errors and in § ^ we introduce 
information content held by an observation with regard to 
the synthesis problem we want to solve. 



3.1 Origin of the errors 

Errors on equivalent widths from spectra of galaxies have 
several origins. The most important one is induced by the 
empirical process of continuum intensity plotting. This leads 
to a wide uncertainty on equivalent widths estimates, caused 
by the fact that the continuum intensity lies on the denom- 
inator in the equivalent width's expression. Contributing to 
the uncertainties in this empirical approach are the blending 
of the absorption lines plus the possible presence of emission 
lines as in active galactic nuclei spectra. The limited wave- 
length range of the study plays also a crucial role in this 
domain. 

Another phenomenon in active galactic nuclei spectra 
is the presence of an additional non stellar continuum which 
leads to the dilution of lines, i.e. the reduction of equivalent 
widths in comparison with those of normal galaxies. The 
effect of this non stellar continuum can be estimated and 
equivalent widths can be corrected; but obviously, this esti- 
mation is also subject to errors. Finally, velocity dispersion 
in the galaxies broadens lines and is also a source of errors. 
Considering these facts, typical errors of 10% are present in 
the equivalent widths of strong lines. 

We did not take into account the errors in the equiva- 
lent width of the lines in the data base itself because they 
are supposed to be negligible compared with errors in the 
equivalent width of t he lin es from the galaxy. However we 
provide the user in § B1.3 with a test able to validate this 



Pi = 



n{w) dw , w — ^p{k) , 



(8) 



hypothesis. If these errors were not negligible, one should 
not perform the synthesis. 



where vr is the a priori probability density function of the 
equivalent widths. We shall evaluate this riA-dimensional in- 
tegral by a Monte-Carlo method (see appendix ^). 

The probability p-, may be zero under two circum- 
stances: (i) Wobs is within the synthetic domain and £^ is re- 
duced to {Wobs} only; (ii) the observation and its associated 
error domain are outside the synthetic domain: nW — 9. 
In each case we have maximum information of different na- 
ture conveyed by the observation: (i) we know that the true 
stellar population lies in the solution set (provided that all 
stars present in the galaxy are also present in the data base); 
(ii) we know that the data base is incomplete, more stars 
must be added. 

In a first approach, it seems reasonable to define the 
information brought by an observation by: I-y = 1 — p^, 
i.e. it is the probability that a stellar population drawn at 
random from the simplex S induces equivalent widths that 
do not fall within the error domain. By 'at random' we mean 
according to n, the a priori distribution of the equivalent 
widths. In practice we use the transform by 93 of an a priori 
distribution of the stellar populations which is uniform over 
S. Defined like this, 1^ is zero if all models are validated 
by the observation (no information) and one if the solution 
set has a zero probability to be reached by a random model 
(maximum information, one cannot get there at random). 

An illustration of the information content within an ob- 
servation is given in Fig|^, see also Fig. 4 of paper II. Here, at 
the 1-a confidence level (7 — 0.683) the information content 
is I-y — 0.92. This observation has therefore a high informa- 
tion content. It is located on a part of the synthetic domain, 
which is seldom reached by a random synthesis following an 
a priori uniform distribution in k. 



4 INFORMATION CONTENT OF AN 
OBSERVATION. 

It is clear that if the synthetic domain W is entirely con- 
tained within the observational error the observation 
cannot discriminate, at the level 7, between the different 
stellar populations. All possible models which fall within £^ 
must be considered as indistinguishable. In other words, the 
observation bears not enough information to differentiate 
between various possible stellar populations of the galaxy. 
The size of the error ellipsoid is in some way the resolu- 
tion at which we try to sort between the different models. 

On the contrary, if the observation is so good that one 
can consider the 'error-bars' to be null, the observation in- 
deed convey a maximum information: Wobs must be syn- 
thesized by the data base. Under this extreme hypothesis, 
£~f is reduced to the single point {Wobs} which has a zero 
measure in W. 

Following a suggestion of R. Barrett we define the in- 
formation brought by the observation as being a function 
of the measure of £^ within W. More precisely will be 
a function of the probability p-y that the image by lys of a 
point drawn at random from the a priori distribution of the 
stellar populations falls within the error ellipsoid. That is: 



5 ERROR ANALYSIS IN PRACTICE. 

In this section, we limit ourselves to the practical results 
of the error analysis. We refer to appendix B for all the 
computations needed to derive the various matrices which 
appear below. 

The key parameter here is the number ris of stars neces- 
sary to define the solution (or the extreme solutions) . In the 
overdetermined case ris = n* , while in the under-determined 
case Tis = riA 4-1 if the galaxy can be synthesized and ris = n\ 
if it is not. It is, in both under- and overdetermined cases, 
less or equal to n\ + 1. We have a regular case if ris = tia -I- 1 
and a singular one if ns < riA -I- 1. We consider separately 
these two cases. In both cases, V stands for the variance- 
covariance matrix of the observations and s for the list of 
the indices of the stars retained in the solution. 

5.1 Regular case ns = riA -I- 1. 

Note that in this case the galaxy belongs to the synthetic 
domain. A solution k is any barycentric combination of uk 
extreme solutions denoted by fc^ (uk may be equal to one). 
At least n* — (nA -1-1) components of ks are equal to zero. The 
components of ks are in addition subject to errors, due to the 
presence of errors in the data. In the tangent approximation 
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Table 1 . The extreme solutions and errors associated with the ob- 
servation and data base illustrated in Fig. |^. Solution number s is 
ks ; uks is the standard deviation vector associated with the solu- 
tion. For the extreme solution Nr.l we also give A = 2ac-y in order 
to allow a direct comparison with Fig. |^. The parameter c~f deter- 
mines the size of the error ellipsoids. We have: Cj = [F~2 i'y)V^'^ ^ 
e.g. with 7 0.683 and nx = 2 degrees of freedom one gets 
Cj ~ 1.515. The quantity A is the difference between the maxi- 
mum and the minimum attained by the stellar contributions of an 
observation which is constrained to move on the error ellipsoid. 
In spite of the non-linear nature of the problem illustrated here, 
it is clear from Fig. ^ that the error analysis under the tangent 
approximation is excellent. 





ki 


k2 




k4 




fcl 


0.306 


0.245 


0.000 


0.449 


0.000 


crfci 


±0.115 


±0.094 




±0.119 




A 


0.347 


0.286 




0.361 




k2 


0.491 


0.000 


0.218 


0.291 


0.000 


ak2 


±0.094 




±0.111 


±0.124 




ks 


0.231 


0.000 


0.000 


0.308 


0.462 


aks 


±0.111 






±0.115 


±0.150 



(small errors) , the variance-covariance matrix of the extreme 
solutions is a singular matrix given by: 



O^O 



K = W 



The elements of the matrix W are: 



= 1, 



w. 



js{i)) 



(9) 



(10) 



where s(i), i — 1, . . . ,ns run over the stars retained in the 
solution. The set of acceptable extreme solutions is an ellip- 
soid centered around the extreme solution ks . Its character- 
istic matrix Ps is given by: 

P. =W^(^r' JJ)w, W = J-^B, 

where J is a + 1 by + 1 diagonal matrix: 

j = 1, . . . ,nx; 
j = nx + l. 



(11) 



The set of acceptable solutions may be considered as the 
convex hull of all these acceptable extreme solutions. 

We have illustrated in Fig. |l| the error analysis per- 
formed on an observation subject to errors of 10%. The pop- 
ulation synthesis was done with the same data base used in 
paper II. The results are given in Table |l|. The fc,, s = 1,2,3 
are the three extreme solutions; aks is the standard devia- 
tion vector associated with the solution (it is the square root 
of the diagonal of V^^); A is the extreme range attained by 
the fciS by a point, which is constrained to move on the er- 
ror ellipsoid. The analysis was done at the confidence level 
of l-cr, that is: 7 ~ 0.683. 



5.1.1 Merit order among extreme solutions 

The product of the non-zero eigenvalues of Vs, is propor- 
tional to the surface of the ellipsoid /C^. This quantity al- 
lows to sort the extreme solutions in order of merit, from 




Figure 1. Illustration of the error analysis for the extreme so- 
lution fcl . Around the observation Wobs = (4, 6) there is an er- 
ror ellipsoid defined by a variance-covariance matrix V. Here this 
matrix is diagonal, the square root of its diagonal elements (the 
standard deviations of the observations) are set to 10% of Wobs- 
The synthetic domain is defined by a data base of 5 stars identi- 
cal to the one given in Table 1 of paper II. This observation and 
its error ellipsoid contain the information Ij 0.92 at the l-cr 
confidence level: 7 = 0.683. 



the smallest surface (highest merit) to the greatest surface 
(smallest merit) . In a synthesis with many extreme solutions 
it is advisable to retain only the solutions which possess the 
highest merit. 



5.2 Singular case < nx + !■ 

In the singular case, an approximate solution is searched by 
minimizing a distance from the observation to the synthetic 
domain. This distance is usually the elliptical distance D, 
defined by Eq. (j^, and it requires a positive definite ma- 
trix S. However, we suppose below that D is the Euclidean 
distance t hrou gh the variable transformation described in 
appendix B2.5. The ellipse has been transformed into the 



unit circle and S = I. 

The error analysis is complicated by the fact that one 
must perform a projection on a plane which is tangent to 
the synthetic domain around the approximate solution. We 
describe step by step the operations which lead to the result. 



Step A: 

ments: [A 



First, construct the nx by Us matrix A of ele- 



Ijs(i)- The index s{i) runs 



over the Us stars defining the solution. We need also the 
diagonal matrix Ja = diag(/syni, • . • , Isynn^^). 

Step Q: The matrix G = J^^A is constructed. A rank de- 
ficiency QR algorithm is performed on G (e.g. using dgeqpf 
from lapack). Following this operation one gets: GII — QR. 
The matrix Q is a by nx square matrix partitioned into 
two blocks: Q = (p|c|). The first block is formed by the ris — l 
first columns of Q. An orthogonal projector H — pp^ is con- 
structed. 
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step V: We define the matrix B equal to A with 
a supplementary line of 'ones' and the matrix K = 
(B B)^^B J where J has already been defined (J — 

diag(/synl, . . .,Lynn^,l) )■ 

The variance-covariance matrix of the solution is: 



K 



HVH' 

0^ 



Step P: We define the oblique projector M: 



M = pp^ + qq^Vp(p^Vp)^p^, 



and compute: W = J ^B. The matrix P entering the defini- 
tion of the acceptance zone in Eq. (fTl) is given by: 



,.,t/M^V^M 



w, w 



Note that in the regular case we have: K — B~^J and in 
the singular case we have: K = (B"^B)~^B"^J. Therefore K is 
always the least square solution of the linear system BK = J. 
(In both cases B is constructed with Wayn and the stars 
contributing to the synthesis.) 



6 THE EXAMPLE OF NGC 3521. 

To illustrate our study we take, as an example, the clus- 
ters of Bica's data base (see Bica 1988) and its galaxy S3 
(NGC 3521). As our aim is to visualize the error zone, we 
shall use only two lines (i.e. n\ — 2) this means two degrees 
of freedom for the error analysis. (Indeed our analysis is valid 
for any number of lines.) From the 35 clusters constituting 
the data base, we select according to Bica's results the eight 
clusters (i.e. n* — 8) contributing to the synthesis. We have 
chosen the two largest equivalent widths of the galaxy which 
correspond to the Call K line {\k ~ 3933A) and the CHG 
band (Ac = 4301 A); the reason of this choice is that the rel- 
ative errors on the strongest lines are generally the smallest. 
Typical errors on such equivalent widths are in the order of 
10%, but for the sake of argument, we take an error of only 
3%. We also assume, that there is no correlation between 
equivalent widths errors (p = 0) even if this is probably not 
the case in reality. 

Fig. 1^ shows that the galaxy can be synthesized ex- 
actl y. O ur method gives 16 different extreme solutions. Ta- 
ble pil shows these several extreme solutions ks where 
s = 1, . . . , 16; also shown are the standard deviations vectors 
crks corresponding to 1-a confidence level; and the standard 
deviation a'ks taking into account a 1% error in the data 
base. The solutions are sorted by increasing order of merit 
(i.e. surface"^). 

Figure |2| demonstrates clearly the impact of the syn- 
thetic domain structure D on the solution errors. Indeed 
the clusters constituting the data base are nearly aligned 
in W-space and their continuum intensities are such that 
the synthetic domain is compressed (i.e the lines joining the 
clusters are very close and nearly parallel to each others). 
As a consequence large errors are present in the extreme 
solutions (see also figure It is important to bear this 
situation in mind when a data base of clusters (or stars) is 
chosen. Ideally the resolution at which the equivalent widths 
space is sampled should be adapted to the quality of the ob- 
servations. Better observations allow a finer coverage of the 



0.0 




Figure 3. The extreme solution of NGC 3521 synthesis that cor- 
responds to the highest merit. The three clusters contributi ng t o 
the solution are number 2, 7 and 15 (see solution fci in table Bl). 
The ellipsoid of errors is the tangent approximation of image by 
application ip^-^ of the error ellipsoid of figure ti. 



VK-space i.e. more clusters (or stars) may be added. We see 



in table 



Bl 



that the contribution of data base errors to the 
solution er rors i s negligible. Finally, the best solution (as 
defined in S 5.1.1) is solution Nr. 1; the solutions 1, 4 and 11 
are the three solutions among the sixteen extreme solutions 
that agree well with Bica's solution (see fcgjj,^^ in Tab. Bl). 
This result is very satisfying, because we used minimal in- 
formation by considering only two equivalent widths. 



7 CONCLUSION. 

We provide explicit formulae to perform a complete error 
analysis on stellar population synthesis by means of a data 
base. The hypothesis underlying the results are that errors 
on the data base itself are negligible compared with errors 
on the observed galaxy. The method also supposes that the 
errors are reasonably Gaussian and further demands the 
knowledge of the error matrix (the variance-covariance ma- 
trix) of the observed equivalent widths. 

The results provided are (i) the variance-covariance ma- 
trix of the extreme solutions (regular case) or of the approxi- 
mate solution (singular case) ; (ii) the matrix defining a zone 
in the solution space where the stellar populations are in- 
distinguishable at a confidence level 7. 

We also introduced information contained in an obser- 
vation with respect to the problem to be solved. This in- 
formation is evaluated by a Monte-Carlo method using as 
input the a priori distribution of the stellar population. We 
further designed a method to generate a uniform a priori 
stellar population distribution. The same method allows us 
to transform the complex hyper-tetrahedron boundary con- 
straints into simple bounds constraints. This transformation 
simplifies considerably the search, when needed, for a mini- 
mum in the stellar population space. 
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Figure 2. The synthetic domain of the 8 Bica's clusters identified by their numbers (2—7,9,15) in the 2 dimensions equivalent widths 
space corresponding to the Call K line and the CH G band. The cross is the position of NGC3521 in this space and the ellipsoid 
corresponds to a 3% standard deviation on thegalactic equivalent widths at a confidence level of 7 = 0.683. At this level the information 
content of this observation is I-y = 0.92 (see §W). 
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APPENDIX A: UNIFORM MAPPING OF A 
SIMPLEX 

In this appendix we show how to transform an hyper-cube 
into an hyper-tetrahedron. This operation is continuous and 
transforms a uniform density on the hyper-cube into a den- 
sity equally uniform on the hyper-tetrahedron. 

The transformation ip allows a simplification of the con- 
straints fe > subject to ^ . fci = 1 (fe G iS) through the 
mapping fc = tp{u) of the variables Ui subject to simple 
bounds constraints of the type: < Ui < 1. At the same 
time the uniform coverage of the hyper-tetrahedron S per- 
mits to evaluate the information integral (^. 

We also want the change of variables tp to be continuous 
and twice derivable, so that 1/) can be used in an optimization 



program where the v aria bles are subject to the constraints 
fc G <S as needed in § 2.3. 

One can think of several ways to map k = 4>{u), but if 
the transformation is not chosen carefully there is a risk that 
7/; selects preferentially certain zones of the k parameters 
space. In order to avoid this possible bias, we would like to 
find a change of variable ensuring that a uniform coverage 
of the u space (the cube) transforms to a uniform coverage 
of the fc space (the tetrahedron). We found two methods 
satisfying these requirements. 



Al Method 1. 

Let us consider the n* independent random variables Ki, 
following all the same exponential distribution. This dis- 
tribution is, for fc > a probability density function: 
/(fc) = exp(— fc). The joint distribution of the multi- 
plet {Ki , Kn^ ) has the density (fci , . . . , fc„ J = 
OI-i /(^O ~ 6xp(— X/T-i ^'^'^ constant on "Y^^l-i = 
a. The conditional density on ^"^^^ fci = a is therefore uni- 
form for any a > 0, in particular for a = 1. 

Now the random variable Ui = ^ f{k) dfc = 1 — 
exp(— isTi) is uniformly distributed on < Mi < 1 so, 
Ki = — ln(l — Ui) or equivalently Ki = —\n{Ui) is expo- 
nentially distributed. In conclusion, the mapping defined 
by: 



- ln(ui 



< < 1 , 



(Al) 



ensures: (i) ki > 0; (ii) a uniform coverage of the ki domain; 
and (iii) that Ui is only subject to simple bounds constraints 
< iti < 1 , i = 1, . . . , n*. 

The only problem with this method is that we use 
parameters while only n*, — 1 are necessary. The search algo- 
rithm may loose time exploring those lines in u space where 



does not change because the ki differ only by a scal- 
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ing factor. We thus designed 'method 2' in order to solve 
this problem. 



A2 Method 2. 

Here we want to map the ki using only n* — 1 parameters 
of type Ui. Again the Ki are exponentially distributed as in 
'method 1'; this ensures, as shown above, that the density 
is uniform on Yl'tli^^ ~ ^■ 

We now introduce the new variables Qi. For the sake 
of argument we limit ourself to a = 4 example, but it 
should be clear that the method works for any n* > 0. We 
define: 

Qi= K4 + K3+K2 + K1, < Qi < 00 ; 
Q2Qi= K4 + K3 + K2, < Q2 < 1 ; 

Q3Q2Qi=Ki + K3, 0<Q3<1; 
QiQ3Q2Qi=K4, 0<Q4<1. 



The Jacobian of this change of variables, J4 



a(fci 



d{qi, ■■■,qA)' 



is given by: J4 = <?i<?293- The four values (Qi, Q2, Qa, Q4) 
have the density: 

f{qi,q2,q3,q4) = e~''^qlqlq3 . (A2) 



This demonstrates that the Qi are independent and have 
densities: 



3! 



hiqi) 

h{q2) 
fsiqs) = 2(73 , 

fiiqi) = 1 , 



3 — Ql 

9ie 



< qi < 00; 
< 92 < 1 ; 
< <?3 < 1 ; 
< 94 < 1 . 



(A3) 



The variable Qi follows a gamma distribution, (as required 
for a sum of independent exponential random variables) and 
the other variables follow a power distribution of decreasing 
index, the last one being uniform. 

Now if Qi is given, say Qi = X/^^i ~ ^' 
maining Q2, Q3 and Q4 cover uniformly (because the Ki 
are exponential) the hyper-tetrahedron defined by (^. We 
have: 



K4 — Q4Q3Q2 , 

K3^{1~ Q4)Q3Q2 , 
K2 = (1-Q3)Q2, 

ATi = 1 - Q2 . 
Finally if we define: 

Un= / nq"^^dq = Q 



(A4) 



(A5) 



the variables [/„ are uniform and independent. We then 

have a power distribution of index n: Qn^-n+i = Un (for 
n = 1, . . . , n* — 1) when Un is uniformly distributed in [0, 1]. 
Therefore the change of variables tp: 



k4 — U1U2 U3 , 
fc3 = (1 - ui)m|m| , 
k2 

fcl = 1 



(1 

(1 -u|)it| , 



(A6) 
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Figure Al. Uniform filling of an equilateral triangle (i.e. an 
hyper-tetrahedron in a 3D space) by means of 'method 2'. The 
points are obtained by using two independent random vari- 
ables C/i and U o un iformly distributed over [0,1]. According to 
Eqs. (A7) and (A8) the change of variables is: fci = 1 — ,Ju2, 



fc2 = (1 - mi)v/m2 and ^3 = ux,/u2 



maps the k space uniformly if ui , U2 and U3 are drawn from 
three independent random variables Ui,U2,U3 uniformly 
distributed in [0, 1]. Now for any > 1, tp is defined by: 



kl — 1 '^J}_-[ ^ , 



n 



and for n — 2, 



, — 1: 



(A7) 



(A8) 



i — ri+ — n-t-1 



where all Ui a re su bject to the constraint: < < 1. We 
present in Fig. Al an example of such a mapping for = 3. 



APPENDIX B: COMPUTATIONS NEEDED 
FOR THE ERROR ANALYSIS. 

We would like to establish a relationship between a small 
variation dW around Wobs and the resulting variation dk 
around the solution k. The variation dk is constrained to 
Er:,rffc,: = Oby E::^k^ = l■ 

Bl The galaxy can be synthesized. 

Wobs belongs to the synthetic domain (Wobs £ W) if the 
galaxy can be synthesized and the many solutions are con- 
tained within the convex hull of a set of extreme solutions 
(see paper II). An extreme solution has by definition at least 
— {nx + 1) components ki equal to zero, what remains is 
the evaluation of -I- 1 components of k. 

Bl.l Error analysis of the extreme solutions. 

We discard the stars with zero contribution and set n* = ns 
(recall that ris = 4- 1). This allows us to keep the same 
notation. Note that the new matrix A is formed of the riA -I- 1 
columns of the original matrix in system (^) corresponding 
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to the stars which have not been set to zero. In the same 
way (p stands here for the original ip restricted to these stars. 

By adding the line to A we form a wa + 1 by riA + 
1 matrix B and exists, because the problem has an 

(extreme) solution. We have: 



fe = B 



B 



A 



(Bl) 



From here we have dBk + Bdk = 0, the matrix dB is dA 
supplemented with a line of zeros. We have dAji = dWjIji 
and [dAfeJj = dWj Iji^i = dWjlsynj- Let us introduce 
the diagonal matrix J, where Jjj — laynj for j — 1, . . . ,nx 
and the last diagonal element is Jn;^+i,n;^+i = 1. We then 
have: 



Bdk 



-dBk 



dW 

J 



(B2) 



As none of the continua Iji are null, exists and we define 
W = J~^B. The matrix W has a very simple form. For ex- 
ample: if we have nx — 2, the number of stars of an extreme 
ray is = 3 and W equals: 



<»'obsl-l-''ll)T 



("'obs 



synl 



^syn2 



(»'obsl-"'l2)T 
("'obs2-"'22)T 



synl 
J22 
syn2 



("'obsl-W'l3)T 



("'obs2-»'23)T 
1 



synl 
J23 



iyn2 



P is usually equal to the inverse of the variance-covariance 
matrix Vfc. The problem is that Vfc, as it is clear from (B4), 
does not possess an inverse. Back to the definition of (p~^{£-y) 
we have: 

IC^ = {k' I Mk') - ^{k)fM-^[^{k') - <^(fc)] < f;/(7)} , 

= {fe' I dW'^M-^dW < F-^\-i)} , 

where we have set dW — ip{k') — ip{k). The above expression 
is equivalent with: 

lC^<-ik'\idW^O){l;'l){'^^)<F-Jm. 

The choice of the extra elements in the central matrix is arbi- 



trary. We chose zero for simplicity. We can now use Eq. (B3) 
and write: 

= {fe' I (fe' - fe)^W^ ( 0^ ' ) -k)< F;^m . 

Therefore: 

,TI V"^ 



p^w.(--)w. 



B1.3 Contributions to errors coming from stellar library 



If we define K = W 
as 



B ^J, equation (B2) can be written 



Wdfc = -( Q j or dfe = -K| 1 



/ 



(B3) 



We can now compute the variance-covariance matrix of the 
fc's. Let Vfc be that matrix. In the tangent approximation 
we have — (dkdk^): 

Vfe = (K(^|j^)(dW^O)K^), 



K 



{dWdW') 0\^T 



0' 







If V designates the variance-covariance matrix of the obser- 
vation we finally have: 



(B4) 



As a bonus, we note that it is possible to express the 
solution fe in terms of K. We have B = JW, B~^ = KJ"^ 



and the solution (Bl) can be expressed by: 



fe = KJ 



Let's call Jwi and J/,; the diagonal matrices such that for 
j = 1, . . . , riA we have [Jwijjj = Ijiki and [Jii]jj = (Wsyn j — 
Wji)ki. In addition we have for j — n\ -t- 1 that [Jw?']jj — 
[Jji]jj ~ 1. Define further Kwi = B~^Jwi and K g = B^^Jn 
Following the same reasoning as in section Bl.l we find that 



the contribution of the stellar library to the total variance- 
covariance matrix of a solution is: 



'H*) 



or 



Wi + Kji 



0^ 



Where \/wi and Vji are the variance-covariance matrices of 
the Wji and Iji of star i. A comparison of \/k(*) and Vfc is 
able to tell if the errors introduced by the stellar library are 
negligible relative to the errors on the observed galaxy. 



B2 The galaxy cannot be synthesized 

An approximate solution Wsyn is found if the galaxy can- 
not be synthesized. A manifold of solutions is defined by a 
number of stars less or equal to the number of equiva- 
lent widths. Here we distinguish ris from the n* number of 
stars in the data base. This manifold is called the synthetic 
'surface' in paper I. 



B1.2 Acceptance region of the extreme solutions 

At the level 7, all points around the solution fe are indistin- 
guishable and give an image close to the observation Wobs. 
The acceptable solutions k' are in a set /C^ called the 'ac- 
ceptance region'. We define: /C^ — (f^^{£-y). 

In the tangent approximation around fe, we have: 

IC-, = {k' I (fe' - fe)^P(fe' - fe) < F-/(7)} , 

where, as in Eq. (^, F^2 is the repartition function of a 
variate possessing nx degrees of freedom. That is also 
the equation of an ellipsoid around fe. The 'weight' matrix 



B2. 1 Projector on the tangent plane. 

If the observational errors dWobs are small around TVobs, 
we can expect small deviations dWay-n around Wsyn. In this 
context, we can approximate the synthetic surface by a plane 



tangent to this surface at Wsyn. The results obtained in Bl 



can be applied to Wsyn since, by definition, Wsyn belongs to 
the synthetic domain. In this context let's A and B stand for 
the matrices A and B, where Wobs have been replaced by 



Wsyn. By (B2) Wsyn + dWsyn belongs to the tangent plane 



if the dk are such that: 
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or dWs, 



(B5) 



with Ja = diag(-^i 



synl ; ■ 



This means that dWg, 



belongs to the linear hull of the columns of J^^A provided 
that '^.dki = 0. We maintain that dWsyn belongs to this 
linear hull even if dki ^ 0. 

In fact ^ . dki = is the equation of an hyperplane 
Eo orthogonal to the subspace of dimension one generated 
by e.g. the vector (1, . . . , 1). The hyperplane So is itself a 
subspace since it includes the origin, therefore any vector dk 
is the sum of a vector dfco G So plus a vector d So. By 
construction the solution fe So since ^ . fci = 1, then dk 
can be written dk = dko + dak, where da is a scalar. Now 
we have: 



J"^B(dfeo + dak) , 
J^^Bdfeo + daJ~^Bfe , 

/ dWsyn \ , f 

[ ) + [da. 



Therefore dWsyn belongs to the linear hull of the columns 
of G = Ja^^ i-^- it belongs to the range of G. 

The range of G is obtained by the rank deficiency QR 
algorithm (see Golub and Van Loan 1996 Chapter 5 §5.4.1). 
According to this algorithm it is possible to write, after a 
permutation 11, GII = QR. The array Q is a by nx 
matrix. The Us — 1 first columns of Q span the range of 
G (i.e. they are an orthonormal basis of the tangent plane) 
and the remaining columns span the subspace orthonormal 
to it. If, following this partition, we write Q — (p|q), one can 
construct the orthogonal projector H on the tangent plane, 
we have: 

H = pp^ . 

The projector H is a by nx matrix of rank Us — 1- 



B2.2 Least square solution of the linearized problem 

One obtains W'^y-^ = Wsyu + dWsyn from W'^^^ = Wobs + 
dWobs by minimizing a 'distance' between these VF^yn and 
Vl^ojjg. This distance is usually defined as a positive definite 
quadratic form, therefore the correction dW sy-n is implicitly 
defined by: 



(W:,bs-VFsy„)^5]-^(W;i 



W' ") = 

' ' svn I 



min{W',^,-W'fll-\W',y,,^W'). (B6) 
W 

Via a change of variables, one can describe the 'w eight' ma- 
trix by the identity matrix (see Sect. B2.5). We here- 
after consider that = I, then the synthesis is obtained 
by the orthogonal projection of Wobs on the synthetic sur- 
face i.e. using H. Consequently one obtains VFgyn by adding 
rfWsyn = H(Woi3s — Wayn) to Wsyii. Taking into account 



that H(Wobs - Ws, 

dWsyn = HdWobs . 



we obtain: 



(B7) 



B2.3 Variance- covariance matrix of k 

The second step is to find a relation between dWsyn and dk 
where dk is restricted to the ris stars which contribute to 



the synthesis. Since, by definition Wsyn £ W, the overde- 
termined system Bfe = (0 l)'^ possesses an exact solution 
in the least-square sense. Therefore k is also the solution 
of the invertible square system B Bfe — B (01)'^. In this 
linear analysis, the perturbed W'^y^ is close to the synthetic 
surface and the normalization constraint b^k = 1 is approx- 

• T ■ ■ T ■ 

imately satisfied. It follows that B dBfc-|-B Bdfe = 0. Again 
dBfe is given by Eq. (B2), where J is the diagonal matrix of 
the synthetic continua J — diag(Jsyni, 
leads to: 



dk : 



-(B b; 



-b"j 



(B8) 



Note that B being of full column rank, the expres- 
sion (B B)~^B is equal to b' ' the Moore-Penrose 
pseudo-inverse of B (for a definition see e.g. Harville 1997 

Chap. 20.1). Let's define here K = (B^B)"^B^J and re- 
calling that dWsyn = H dWobs we obtain the variance- 
covariance matrix of the synthesis: 



(B9) 



where V = {dWobsdW'^^^^) is the variance-covariance matrix 
of the observations. 



B2.4 Acceptance region of the least-square solutions 

An alternative solution fe' is considered acceptable, at the 
level 7, if there exists a Wo^s within the 'error' region £^ 
around Wobs such that HWobs ~ Wsyn ~ f{k')- Accord- 
ing to this definition, fe' is acceptable if it belongs to the 
following set: 



JC-, = {fe' I dWiTbsV-^Wobs < f 0^(7)} , 



HdWoi 



dWe, 



(BIO) 



where V is the variance-covariance matrix of the observa- 
tional errors. 

Among all the solutions of HdWobs ~ dWsyn we choose 
the one which is 'closest' to Wobs in the V~^-norm. This 
ensures that the solution is within The minimum of 
dWobsV~^dWobs, which satisfies the constraint HdWobs ~ 
dWsyn is given by: 

dWobs = M dWsyn , 

M = pp^-f qq^Vp(p^Vp) ip^. 

The matrix M is an oblique projector on the regression plane 
defined by the constraint. We indeed have: 



/ dWobs 
V 



MdWs, 




(BU) 



V 



One carmot invert Eq. (BJ ) in order to introduce dk, but 
Eq. (|B5|) sine 

=yn ^= _J-iBdfc. 



we can use Eq. (B5) since K.^ is in the set where dki ~ 
/dW 



(B12) 



Sect. § (|B1.2|), we find 



Combin ing E qs. (B12) and (Bll) and proceeding like in 

(B13) 



{fe'|(fc'-fe)^P(fe'-fc) <F-/(7)}, 
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0^ 









W, W = J"'B. 



(B14) 



B2.5 Weighted least square 

We consider now the case where the distance between two 
points of the W-space is not the Euclidean distance but 
an elliptical distance. This elliptical distance is subordi- 
nated to a symmetric definite positive matrix S. We have 
ds{Wi,W2) = {Wi - W2f'E'\Wi - W2). Usually one 
chooses S equal to V the variance-covariance matrix of the 
data, but it may not be always the case. 

A substitution of variables transforms the elliptical dis- 
tance into a spherical (i.e. Euclidean) one. We want in other 
words that d{Wi,W2) ~ dE(VKio, W2o),where d is the 
Euclidean distance. If Wo stands for the original variables 
and W for the new variables, we define the matrix N by: 
Wo = NW. It now follows that 



d{Wi,W2) = (Wi - W: 



^N{Wi - W2) ■ 



This imposes the condition N"^S~^N = I, which is satisfied if 
we choose N as the matrix appearing in the Choleski factor- 
ization of S (see Golub & Van Loan 1996, section 4.2.1). We 
have S = NN"^, where N is lower triangular. It is straight- 
forward to verify that N satisfies the condition. 

Finally we need to compute the varia nce-covariance ma- 
trix V appearing in Eqs. ( |BE| ) and ( B14 ) from the variance- 
covariance matrix Vo of the original variables. We have 



V = {dWdW^} = 



N-^{dWoWo >N~i^, therefore: 



(B15) 
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Table Bl. The Bica solution and our 16 extreme solutions of galaxy NGC3521 synthesis using as data base the 8 clusters (2-7,9,15) 
retain by E. Bica in his synthesis. Under column fej is the contribution to the galactic luminosity of cluster Nr. i; crks, (s = 1, 2, ...16) is the 
standard deviation corresponding to the extreme solution ks supposing only errors on galactic data ajid cr'ks is the standard deviation 
corresponding to the same extreme solution taking into a<;count errors on data base. Finally, 'surface' is the value of the ellipsoid's surface 
for each extreme solution. The solutions are sorted a<;cording to a decreasing order of merit (increasing order of surfa<;e). 

k2 ks ki fcs fee kr kg feis surface 

0.070 0.060 0.050 0.080 0.060 

0.000 0.000 0.091 0.000 0.348 6.46210-5 

±0.407 

±0.418 

0.000 0.000 0.075 0.000 0.050 7.779 lO'^ 

±0.498 
±0.520 

0.000 0.000 0.009 0.000 0.000 8.75910-5 
±0.132 
±0.136 

0.000 0.146 0.000 0.000 0.284 1.45610"* 
±0.104 ±0.438 
±0.106 ±0.448 

0.086 0.000 0.046 0.000 0.000 2.08510-* 
±0.842 
±0.880 

0.000 0.000 0.165 0.568 0.000 2.12110-* 
±0.032 ±0.700 
±0.033 ±0.723 

0.000 0.105 0.010 0.000 0.000 3.09810" 
±1.029 
±1.077 

0.000 0.231 0.000 0.403 0.000 3.16210"* 
±0.042 
±0.043 

0.000 0.017 0.000 0.000 0.000 3.26110"* 
±0.255 
±0.262 

0.000 0.000 0.061 0.000 0.000 3.76910"* 
±0.190 
±0.198 

0.268 0.000 0.000 0.000 0.186 3.98610"* 
±0.180 ±0.480 
±0.183 ±0.490 

0.352 0.000 0.000 0.220 0.000 5.45810"* 
±0.060 
±0.061 

0.000 0.122 0.000 0.000 0.000 6.74610"* 
±0.117 
±0.122 

0.000 0.000 0.097 0.147 0.000 8.08710"* 
±0.178 ±1.538 
±0.183 ±1.589 

0.044 0.000 0.000 0.000 0.000 2.269 10"^ 
±0.671 
±0.689 

0.236 0.000 0.000 0.000 0.000 2.604 10"^ 

±0.228 
±0.2:-! 1 
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